home *** CD-ROM | disk | FTP | other *** search
/ Space & Astronomy / Space and Astronomy (October 1993).iso / mac / programs / satrack / SS4EXE.ZIP / SEESAT4.ZIP / ASTRO.C next >
Text File  |  1992-01-05  |  17KB  |  632 lines
  1. /*
  2. ASTRO.C
  3. by Paul S. Hirose, 1991 Oct 19
  4. Coordinate transformation and astronomical functions for SEESAT satellite
  5. tracking program.
  6.  
  7. This file is in the public domain.
  8.  
  9. 200 - 299
  10. */
  11.  
  12. #include "B:SEESAT.H"    /* global header */
  13.  
  14. /* library functions:
  15.     asin atan atan2 atof cos log10 printf sin sqrt tan
  16. are used in this file. */
  17.  
  18.  
  19. #if ECOC
  20. extern void printf();
  21. extern double asin(), atan(), atan2(), atof(), cos(),
  22.     log10(), sin(), sqrt(), tan();
  23. #endif
  24.  
  25. #if ECOC | LASERC
  26. extern double dint();
  27. #endif
  28.  
  29. #if STDC
  30. static double dint(x)
  31. double x;
  32. {
  33.     if (x >= 0.)
  34.         return floor(x);
  35.     else
  36.         return ceil(x);
  37. }
  38. #endif
  39.  
  40. #define FLAT 3.35278e-3        /* flattening of earth.
  41.                 WGS-72 value = 3.35278e-3 */
  42.  
  43. /* REFEP is the default epoch to which R.A./dec. will be precessed.
  44. E.g. (per Meeus p. 101):
  45. 3477629251. = epoch 1900.0, 3503926689. = 1950.0, 3530224800. = 2000.0.
  46. */
  47.  
  48. #if ENPRE        /* precession enabled */
  49. #define REFEP 3530224800.
  50. #endif
  51.  
  52.  
  53. /* Bias to equalize SEESAT's internal magnitude scale with the absolute
  54. satellite magnitudes being supplied from the element file or manually.  If
  55. the absolute mags are for a satellite at R kilometers and P percent
  56. illuminated, set MAGOFF to 2.5 * log10((6378/R)^2 * (P/50)).  E.g., for
  57. Molczan's standard conditions of 1000 km and 50%, MAGOFF = 4.02. */
  58.  
  59. #define MAGOFF 4.02
  60.  
  61.  
  62. /* 0 for normal operation.  Non-zero will make static functions and data
  63. extern */
  64. #define TEST 0
  65.  
  66.  
  67. /*######################## LOCAL FUNCTIONS ########################*/
  68.  
  69. #if TEST
  70. #define STATIC        /* precedes static func definitions */
  71. #define SC extern    /* precedes static func declarations */
  72.  
  73. #else
  74. #define STATIC static
  75. #define SC static
  76. #endif
  77.  
  78. SC void eqhor();        /* convert rectangular to az/el, print */
  79. SC void eceq();            /* ecliptical to equatorial coordinates */
  80. SC void point();        /* align +z axis with a vector */
  81. SC void recsph();        /* rectangular to spherical coordinates */
  82. SC void sun();            /* xyz of sun (interpol.) */
  83. SC double sunlon();        /* mean longitude of sun */
  84. SC void sunref();        /* xyz of sun (non-interpol.) */
  85. SC void topos();        /* rectangular coordinates of observer */
  86. SC void yrot();            /* y-rotate vector */
  87. SC void zrot();            /* z-rotate vector */
  88.  
  89. /*########################### LOCAL DATA ###########################*/
  90.  
  91. #if ENPRE
  92.  
  93. /* dircos[] is a rotation matrix.  Initialized by inpre() & used to precess
  94. satellite Right Ascension & declination.  terep is the epoch to which the
  95. coordinates will be precessed. */
  96.  
  97. STATIC struct {
  98.     double ix, iy, iz, jx, jy, jz, kx, ky, kz;
  99. } dircos;
  100. STATIC double terep = REFEP;
  101.  
  102. #endif
  103.  
  104. /* data for observer's location */
  105. STATIC struct {
  106.     double
  107.     h,        /* altitude above sea level */
  108.     lambda,        /* longitude */
  109.     sinlat, coslat,    /* latitude sin & cos */
  110.     zg,        /* north distance, observer to equatorial plane */
  111.     xc;        /* distance, observer to polar axis */
  112. } obs;
  113.  
  114. /* Sin & cos of epsilon, the obliquity of ecliptic.  Since epsilon changes
  115. but .013 deg/century, it's sufficient to use a fixed value, that of 1975
  116. (23.442 deg).  That will strike a balance between 1950 & 2000 epochs. */
  117. STATIC double sineps = .39783;
  118. STATIC double coseps = .91746;
  119.  
  120. /*############################## CODE ##############################*/
  121.  
  122. STATIC void
  123. eqhor(equa, t)
  124. struct vector *equa;    /* vector to object */
  125. double t;        /* epoch */
  126. /* Converts vector *equa from a rectangular (Aries, equator, north) system to
  127. a polar (azimuth, elevation, range) system at epoch t.  Prints the azimuth &
  128. elevation rounded to nearest deg.  No correction for refraction or geocentric
  129. parallax. */
  130. {
  131.     double    lhaa;    /* local hour angle of Aries */
  132.  
  133.     /* rotate to a south, east, zenith system */
  134.     lhaa = thetag(t) + obs.lambda;
  135.     zrot(equa, sin(lhaa), cos(lhaa));
  136.     yrot(equa, obs.coslat, obs.sinlat);
  137.  
  138.     equa->x = -equa->x;    /* makes az run correct direction */
  139.     recsph(equa);        /* convert to spherical coordinates */
  140.  
  141.     printf("az=%d ", (int) (equa->x * ra2de + .5));
  142.     printf("el=%d\n", (int) (equa->y * ra2de +
  143.       (equa->y >= 0. ? .5 : -.5)));
  144. }
  145.  
  146.  
  147. void
  148. dusk()
  149. /* Print azimuth and elevation of the sun.  The next argument(s) on
  150. the command line are taken as the epoch. */
  151. {
  152.     double t;
  153.     struct vector csun;    /* equatorial coordinates of sun */
  154.  
  155.     t = tokmin();        /* get epoch from cmd line */
  156.     sun(&csun, t);        /* get sun's position */
  157.     eqhor(&csun, t);    /* convert to az/el, print */
  158. }
  159.  
  160.  
  161. STATIC void
  162. eceq(eclip)
  163. struct vector *eclip;    /* ecliptical coordinates */
  164. /* Coverts *eclip to equatorial coordinates, i.e., x-rotates by the negative
  165. obliquity of the ecliptic. */
  166. {
  167.     double tempy;
  168.  
  169.     tempy    = eclip->y * coseps - eclip->z * sineps;
  170.     eclip->z = eclip->y * sineps + eclip->z * coseps;
  171.     eclip->y = tempy;
  172. }
  173.  
  174.  
  175. double
  176. fmod2p(x)
  177. double x;
  178. /* Reduces x to range 0 - 2pi */
  179. {
  180.     x /= twopi;
  181.     x = (x - dint(x)) * twopi;
  182.     if (x < 0.)
  183.         return x + twopi;
  184.     else
  185.         return x;
  186. }
  187.  
  188.  
  189. #if ENPRE    /* precession code enabled */
  190.  
  191. void
  192. inpre()
  193. /* Initialize rotation matrix dircos to precess from "epoch" (a global
  194. variable representing epoch of the satellite elements) to terep.  terep is an
  195. external variable in this file.  Assumes that the effect of precession is a
  196. uniform 50".4 increase of longitude per year. */
  197. {
  198.     double a, sina, cosa;
  199.     struct vector *vecp;
  200.     int i;
  201.  
  202.     a = (terep - epoch) * -4.65e-10;    /* luni-solar precession */
  203.     sina = sin(a);
  204.     cosa = cos(a);
  205.  
  206.     /* Ecliptical components of (I, J, K) unit vectors */
  207.     dircos.ix = 1.;        dircos.iy = 0.;        dircos.iz = 0.;
  208.     dircos.jx = 0.;        dircos.jy = coseps;    dircos.jz = -sineps;
  209.     dircos.kx = 0.;        dircos.ky = sineps;    dircos.kz = coseps;
  210.  
  211.     /* Rotate each vector in longitude by the amount of luni-solar
  212.     precession.  Then convert from ecliptical to equatorial coordinates.
  213.     We end up with the old unit vectors expressed in terms of the new
  214.     (i.e., precessed) equatorial coordinate system. */
  215.  
  216.     for (vecp = (struct vector *) &dircos, i = 3; i--; ++vecp) {
  217.         zrot(vecp, sina, cosa);
  218.         eceq(vecp);
  219. }    }
  220.  
  221. #endif
  222.  
  223.  
  224. void
  225. moon()
  226. /* Prints the azimuth, elevation, and % of illumination of the moon.  Epoch
  227. is taken from the next arguments on the command line.  Position is good
  228. within a couple degrees or so. */
  229. {
  230.     struct vector luna, sol, terra;
  231.     double t, t1900, lp, mp, f, lamb, beta, cosb;
  232.  
  233.     t = tokmin();            /* obtain epoch from cmd line */
  234.  
  235.     /* formulas from Meeus */
  236.     t1900 = t - 3477628800.;        /* 1900 Jan 0 12h UT */
  237.     lp = fmod2p(4.72 + 1.5970243e-4 * t1900);    /* longitude */
  238.     mp = fmod2p(5.168 + 1.5835217e-4 * t1900);    /* mean anomaly */
  239.     f = fmod2p(.1964 + 1.6034425e-4 * t1900);    /* dist fm asc node */
  240.     lamb = lp + .1097 * sin(mp);            /* longitude */
  241.     beta = .0895 * sin(f);                /* latitude */
  242.  
  243.     /* get ecliptical components of unit vector to moon */
  244.     luna.z = sin(beta);
  245.     cosb = cos(beta);
  246.     luna.y = sin(lamb) * cosb;
  247.     luna.x = cos(lamb) * cosb;
  248.  
  249.     eceq(&luna);        /* ecliptical to equatorial */
  250.  
  251.     /* terra = unit vector moon to earth */
  252.     terra.x = -luna.x;
  253.     terra.y = -luna.y;
  254.     terra.z = -luna.z;
  255.  
  256.     eqhor(&luna, t);    /* print az el */
  257.  
  258.     /* Get moon->sun unit vector.  sun(), which yields the earth->sun
  259.     vector, comes close enough.   Rotate axes of terra to point +z axis
  260.     at sun.  This will make terra.z equal to the cos of the moon's phase
  261.     angle.  Plug that into Meeus' illumination formula, print. */
  262.  
  263.     sun(&sol, t);
  264.     point(&terra, &sol);
  265.  
  266.     printf("%d%% illuminated\n", (int) (50. * (1. + terra.z)));
  267. }
  268.  
  269.  
  270. void
  271. parall()
  272. /* Prints the parallactic angle:  the angle (at the satellite) between
  273. celestial north and observer's zenith.  The next argument(s) on the command
  274. line are taken as the desired epoch. */
  275. {
  276.     struct vector zenith, topsat;
  277.     double lhaa, t, coslha, sinlha;
  278.  
  279.     t = tokmin();
  280.     lhaa = thetag(t) + obs.lambda;    /* local hr angle of Aries */
  281.     sinlha = sin(lhaa);
  282.     coslha = cos(lhaa);
  283.  
  284.     /* Get equatorial coordinates of satellite at time t. */
  285.     MODEL(t - epoch - toffs);
  286.  
  287.     /* get topocentric sat coordinates */
  288.     topsat.x = sat.x - obs.xc * coslha;
  289.     topsat.y = sat.y - obs.xc * sinlha;
  290.     topsat.z = sat.z + obs.zg;
  291.  
  292.     /* Load struct "zenith" with a unit vector directed from the observer
  293.     to the zenith, expressed in equatorial coordinates. */
  294.  
  295.     zenith.z = obs.sinlat;
  296.     zenith.x = coslha * obs.coslat;
  297.     zenith.y = sinlha * obs.coslat;
  298.  
  299.     /* z-rotate the axes of "zenith" so satellite is in the x-z plane of
  300.     "zenith" (+x half plane).  Then y-rotate axes to point +z axis at
  301.     satellite. */
  302.  
  303.     point(&zenith, &topsat);
  304.  
  305.     /* The +z axis passes through the satellite.  The north pole
  306.     has a y-coordinate of 0 and a negative x-coordinate.  The
  307.     origin of the coordinate system is still at the observer. */
  308.  
  309.     printf("parallactic angle = %d\n",
  310.       (int) (fmod2p(atan2(zenith.y, -zenith.x)) * ra2de));
  311. }
  312.  
  313.  
  314. STATIC void
  315. point(v, d)
  316. struct vector *v, *d;
  317. /* rotates the axes of vector v to align the +z axis with vector d. */
  318. {
  319.     double fz, fz2, r;
  320.  
  321.     fz2 = d->x * d->x + d->y * d->y;
  322.     fz = sqrt(fz2);
  323.     zrot(v, d->y / fz, d->x / fz);
  324.     r = sqrt(fz2 + d->z * d->z);
  325.     yrot(v, fz / r, d->z / r);
  326. }
  327.  
  328.  
  329. STATIC void
  330. recsph(v)
  331. struct vector *v;
  332. /* Converts vector v from rectangular to spherical coordinates.  On exit,
  333. v->x = angle in the x-y plane, v->y = angle with respect to the x-y plane,
  334. and v->z = radius.  Angle x is 0 to 2pi, y is -pi/2 to pi/2. */
  335. {
  336.     double temp;
  337.  
  338.     temp = sqrt(v->x * v->x + v->y * v->y + v->z * v->z);
  339.     v->x = atan2(v->y, v->x);
  340.     if (v->x < 0.)
  341.         v->x += twopi;
  342.     v->y = asin(v->z / temp);
  343.     v->z = temp;
  344. }
  345.  
  346.  
  347. #if ENPRE    /* precession code enabled */
  348.  
  349. void
  350. setep()
  351. /* Sets rotation matrix "dircos" to precess from "epoch" (a global variable
  352. representing epoch of the satellite orbital elements) to the date/time given
  353. on the command line. */
  354. {
  355.     terep = tokmin();    /* get epoch from cmd line */
  356.     inpre();    /* reinitialize rotation matrix */
  357. }
  358.  
  359. #endif
  360.  
  361.  
  362. void
  363. seth()
  364. /* input observer's height in km, recompute rectangular coordinates */
  365. {
  366.     obs.h = atof(*tokp) / xkmper;
  367.     topos();
  368. }
  369.  
  370.  
  371. void
  372. setlat()
  373. /* input observer's latitude, recompute rectangular coordinates */
  374. {
  375.     obs.coslat = atof(*tokp) * de2ra;
  376.     obs.sinlat = sin(obs.coslat);
  377.     obs.coslat = cos(obs.coslat);
  378.     topos();
  379. }
  380.  
  381.  
  382. void
  383. setlon()
  384. /* input observer's longitude, recompute rectangular coordinates */
  385. {
  386.     obs.lambda = atof(*tokp) * de2ra;
  387.     topos();
  388. }
  389.  
  390.  
  391. STATIC void
  392. sun(pos, ep)
  393. struct vector *pos;
  394. double ep;
  395. {
  396. /* Equatorial xyz components of a unit vector to sun at epoch "ep" are
  397. returned in struct "pos".  Accuracy about .1 deg. */
  398.  
  399.     static double t0,    /* start epoch of a 6-day interval */
  400.         dx, dy, dz;    /* change in x, y, z during 6 days */
  401.     static struct vector sun0;    /* sun position at t0 */
  402.  
  403.     double t;    /* e.g. .5 if ep = t0 + 3 days */
  404.  
  405.     if ((t = (ep - t0) / 8640.) > 1.0  ||  t < 0.) {
  406.     /* Position was requested for a time outside the current 6-day
  407.     interval.  So set up a new 6-day interval centered at "ep". */
  408.         struct vector sun1;
  409.  
  410.         t0 = ep - 4320.;    /* 3 days before "ep" */
  411.         sunref(&sun0, t0);
  412.         sunref(&sun1, ep + 4320.);    /* 3 days after "ep" */
  413.         dx = sun1.x - sun0.x;
  414.         dy = sun1.y - sun0.y;
  415.         dz = sun1.z - sun0.z;
  416.         t = .5;
  417.     }
  418.     /* get position by linear interpolation in the 6-day period */
  419.     pos->x = sun0.x + t * dx;
  420.     pos->y = sun0.y + t * dy;
  421.     pos->z = sun0.z + t * dz;
  422. }
  423.  
  424.  
  425. STATIC double
  426. sunlon(epoch)
  427. double epoch;
  428. /* Return geometric mean longitude of sun at "epoch", accurate to about .1
  429. deg.  Formulas adapted from Meeus.  The returned value may be off by a small
  430. multiple of 2pi.  That should cause no problem since in this program the
  431. longitude is only used by sin() & cos(). */
  432. {
  433.     static double e = .016709;    /* eccentricity year 2000 */
  434.  
  435.     double t, l, m, enew, eold, v;
  436.  
  437.     t = epoch - 3477628800.;    /* since 1900 Jan 0 12h */
  438.  
  439.     /* geometric mean long. */
  440.     l = fmod2p(4.881628 + t * 1.1946383e-5);
  441.  
  442.     /* mean anomaly */
  443.     m = fmod2p(6.256584 + t * 1.1945812e-5);
  444.  
  445.     /* solve Kepler's equation to .1 deg. */
  446.     enew = m;
  447.     do {
  448.         eold = enew;
  449.         enew = m + e * sin(eold);
  450.     } while (FABS(enew - eold) >= 1.7e-3);
  451.  
  452.     /* tan is asymptotic at pi/2, so if eccentric anomaly is within .1
  453.     deg of pi radians, we'll just make true anomaly = pi. */
  454.     if (FABS(enew - pi) <= 1.7e-3)
  455.         v = pi;
  456.     else
  457.         v = 2. * atan(sqrt((1. + e) / (1. - e)) * tan(enew / 2.));
  458.  
  459.     return l + v - m;
  460. }
  461.  
  462.  
  463. STATIC void
  464. sunref(pos, ep)
  465. struct vector *pos;
  466. double ep;
  467. /* Put the mean equatorial xyz components of a unit vector to sun at
  468. epoch "ep" into struct "pos".  Good to about .1 deg. */
  469. {
  470.     double lon, sinlon;
  471.  
  472.     lon = sunlon(ep);
  473.     sinlon = sin(lon);
  474.  
  475.     pos->x = cos(lon);
  476.     pos->y = sinlon * coseps;
  477.     pos->z = sinlon * sineps;
  478. }
  479.  
  480.  
  481. double
  482. thetag(ep)
  483. double ep;    /* epoch */
  484. /* Returns Greenwich hour angle of the mean equinox at ep.  If ep is UT1,
  485. returned value is within .001 sec of time compared to the formula in the '85
  486. Astronomical Almanac, for 1990 through 2010.  As a side effect, this function
  487. sets ds50 to days since 1950 Jan 0 0h UT. */
  488. {
  489.     double theta;
  490.  
  491.     ds50 = (ep - 3503925360.) / xmnpda;
  492.     theta = .27524987 + 1.00273790935 * ds50;    /* revolutions */
  493.  
  494.     return (theta - dint(theta)) * twopi;
  495. }
  496.  
  497.  
  498. STATIC void
  499. topos()
  500. /* Compute rectangular coordinates of observer on the ellipsoid.
  501. Formulas from Baker & Makemson. */
  502. {
  503.     double c, s;
  504.     static double f = FLAT;        /* flattening of earth */
  505.  
  506.     c = 1. / sqrt(1. - (f + f - f * f) * obs.sinlat * obs.sinlat);
  507.     s = c * (1. - f) * (1. - f);
  508.     obs.xc = (c + obs.h) * obs.coslat;
  509.     obs.zg = -(s + obs.h) * obs.sinlat;
  510. }
  511.  
  512.  
  513. int
  514. xyztop(t)
  515. double t;    /* epoch */
  516. /* If satellite is below horizon and aflag = 0, xyztop() returns 0.
  517. If satellite is above horizon or aflag is true, loads globals
  518. elsusa, azel, radec and latlon with appropriate values and returns 1. */
  519. {
  520.     struct vector suneq;    /* equatorial sun position */
  521.  
  522.     double temp, tempx, tempy, lhaa, sinlha, coslha;
  523.  
  524.     lhaa = thetag(t) + obs.lambda;        /* local hr angle Aries */
  525.     sinlha = sin(lhaa);
  526.     coslha = cos(lhaa);
  527.  
  528.     /* load azel & radec with sat coordinates translated to observer */
  529.     radec.x = azel.x = sat.x - obs.xc * coslha;
  530.     radec.y = azel.y = sat.y - obs.xc * sinlha;
  531.     radec.z = azel.z = sat.z + obs.zg;
  532.  
  533.     /* Rotate azel into a south/east/zenith system */
  534.     zrot(&azel, sinlha, coslha);
  535.     yrot(&azel, obs.coslat, obs.sinlat);
  536.  
  537.     if (azel.z < 0. && aflag == 0)
  538.         return 0;    /* below horizon and "all" flag not set */
  539.  
  540.     azel.x = -azel.x;    /* makes az go correct direction */
  541.     FTEST((200, 3, 4, &azel.x, &azel.y, &azel.z));
  542.     recsph(&azel);        /* convert to spherical */
  543.  
  544. #if ENPRE    /* precess RA & dec */
  545.     tempx   = radec.x * dircos.ix + radec.y * dircos.jx +
  546.       radec.z * dircos.kx;
  547.     tempy   = radec.x * dircos.iy + radec.y * dircos.jy +
  548.       radec.z * dircos.ky;
  549.     radec.z = radec.x * dircos.iz + radec.y * dircos.jz +
  550.       radec.z * dircos.kz;
  551.  
  552.     radec.x = tempx;
  553.     radec.y = tempy;
  554. #endif
  555.     FTEST((201, 3, 4, &radec.x, &radec.y, &radec.z));
  556.  
  557.     sun(&suneq, t);         /* suneq = xyz of sun */
  558.  
  559.     /* get dot product of observer->sat and observer->sun vectors */
  560.     temp = radec.x * suneq.x + radec.y * suneq.y + radec.z * suneq.z;
  561.  
  562.     recsph(&radec);        /* convert sat xyz to RA, dec, range */
  563.  
  564.     /* Illuminated fraction of sat = (1 - cos(i)) / 2, where i is phase
  565.     (i.e., illum. angle) of sat.  Apparent magnitude is a function of
  566.     illuminated fraction, range squared, absolute magnitude, and a scale
  567.     offset factor.  We can simplify considerably by observing that:  1)
  568.     cos(i) = temp / (|sat| * |sun|) (this is a restatement of the
  569.     definition of a dot product), 2) |sat| (range of satellite) is now in
  570.     radec.z, 3) |sun| = 1 since sun() always generates a unit vector.  In
  571.     addition, we can drop the "/ 2", adjusting MAGOFF accordingly. */
  572.  
  573.     apmag = 2.5 * log10(radec.z * radec.z / (1. - temp / radec.z))
  574.       + abmag + MAGOFF;
  575.  
  576.     /* latitude & longitude */
  577.     latlon.x = sat.x;
  578.     latlon.y = sat.y;
  579.     latlon.z = sat.z;
  580.     recsph(&latlon);    /* convert to spherical coordinates */
  581.  
  582.     /* figure longitude relative to local meridian or Greenwich,
  583.     depending on greenwich flag. */
  584.     latlon.x = fmod2p(latlon.x - lhaa + (gflag ? obs.lambda : 0.));
  585.  
  586.     if (latlon.x > pi)
  587.         latlon.x -= twopi;
  588.  
  589.     latlon.z -= mean_r;    /* mean_r = mean radius of earth */
  590.  
  591.     /* Get sun elevation at satellite by rotating sun coordinate axes
  592.     so the satellite is on the positive z-axis.  Sun elevation = arc
  593.     sin z.  Add dip of horizon due to height of satellite. */
  594.  
  595.     point(&suneq, &sat);
  596.  
  597.     temp = (asin(suneq.z) + atan(sqrt((2. + latlon.z) * latlon.z)))
  598.       * ra2de;
  599.     FTEST((202, 1, 1, &temp));
  600.     elsusa = temp + ((temp >= 0.) ? .5 : -.5);
  601.  
  602.     return 1;    /* signifies sat elev >= 0 */
  603. }
  604.  
  605.  
  606. STATIC void
  607. yrot(vec, sint, cost)
  608. struct vector *vec;
  609. double sint, cost;
  610. /* y-rotates coordinate axes for "vec" by angle t, where "sint"
  611. and "cost" are the sine and cosine of t. */
  612. {
  613.     double tempz;
  614.     tempz  = vec->x * sint + vec->z * cost;
  615.     vec->x = vec->x * cost - vec->z * sint;
  616.     vec->z = tempz;
  617. }
  618.  
  619.  
  620. STATIC void
  621. zrot(vec, sint, cost)
  622. struct vector *vec;
  623. double sint, cost;
  624. /* similar to yrot(), except rotate about z axis */
  625. {
  626.     double tempx;
  627.  
  628.     tempx  = vec->y * sint + vec->x * cost;
  629.     vec->y = vec->y * cost - vec->x * sint;
  630.     vec->x = tempx;
  631. }
  632. ean